/******************************************************************************
 *	COPYRIGHT (C), 1999-2000 SEIKO EPSON CORP.
 *  ALL RIGHTS RESERVED
 *
 *	File: ATAN2.C
 *
 *	C Standard Library Function "atan2"
 *	The "atan2" function calculates the arc tangent of "y/x".
 *
 *	Revision History
 *	12/08/1999	Created by H.Ishida
 *  01/06/2000  Fixed a bug by H.Ishida
 *  02/14/2000	change value M.Igarashi
 *
 * -- Format of floating point (double precision) ---
 * 
 * 63|62         52|51               32|31                                0
 *  -----------------------------------------------------------------------
 * | |   exponent  |                    fraction                           |
 *  -----------------------------------------------------------------------
 * |    12 bits          20 bits       |             32 bits               |
 * |            lower word             |            higher word            |
 ******************************************************************************/

#include <math.h>
#include <smcvals.h>

#define	K0		0.445413400592906803197511e2
#define	K1		0.77477687719204208616481e2
#define	K2		0.40969264832102256374186e2
#define	K3		0.666057901700926265753e1
#define	K4		0.1589740288482307048e0
#define	S0		0.445413400592906804445995e2
#define	S1		0.92324801072300974840693e2
#define	S2		0.62835930511032376833267e2
#define	S3		0.1550397755142198752523e2

double _atan(double );

/******************************************************************************
 *	Function: double atan2(double dfX, double dfY)
 ******************************************************************************/
double atan2(double dfX, double dfY)
{
	long lX,lX2;
	unsigned long ulLx,ulHx,ulHx2;
	double	dfTemp;
	double dfX2,dfRet,dfX3;
	
	
	GETW_L(lX,dfX);		// get low
	GETW_L(lX2,dfY);		// get low
	GETW_H(ulHx,dfX);		// get high
	GETW_H(ulHx2,dfY);		// get high
	
	SET_SIGN( dfX2, dfX );	// clear sign

	
	dfTemp=dfX+dfY;
	if (dfTemp==dfX){
//	if((dfX+dfY)==dfX){
		dfRet=HLFP;
		if( (lX&0x80000000) == 0 || ( lX ==0x0 && ulHx==0x0 )){
			return dfRet;
		} else {
			return -dfRet;
		}
	}else if( (lX2&0x80000000) != 0 || ((lX2&0x80000000) != 0 ) && (ulHx2 != 0)) {
		dfRet = WHP-_atan(dfX2/dfY);
		if( (lX&0x80000000) == 0 || ( lX ==0x0 && ulHx==0x0 )){
			return dfRet;
		} else {
			return -dfRet;
		}
	}
	else if((lX&0x80000000) == 0 || ((lX&0x80000000) == 0) && (ulHx !=0x0 )){
		dfRet = _atan(dfX/dfY);
		return dfRet;
	} else {
		dfRet = -_atan(-dfX/dfY);
		return dfRet;
	}
}

/******************************************************************************
 *	Function: double _atan(double dfX)
 ******************************************************************************/
double _atan(double dfX)
{
	long lX;
	unsigned long ulHx;
	double dfXsq,dfRet,dfTemp;
	
	GETW_L(lX,dfX);		// get low
	GETW_H(ulHx,dfX);		// get high
	
	if( (lX&0x80000000) != 0 || lX < 0x3fda8279 || (lX == 0x3fda8279 && ulHx < 0x99fcef32 )){	// check qarea
		dfXsq = dfX*dfX;
		dfRet=K4*dfXsq+K3;
		dfRet=dfRet*dfXsq+K2;
		dfRet=dfRet*dfXsq+K1;
		dfRet=dfRet*dfXsq+K0;
		dfTemp=dfXsq+S3;
		dfTemp=dfTemp*dfXsq+S2;
		dfTemp=dfTemp*dfXsq+S1;
		dfTemp=dfTemp*dfXsq+S0;
		dfRet = dfRet/dfTemp*dfX;
		return dfRet;
	}else if( ((lX&0x80000000) == 0 )&& ( lX > 0x4003504f || (lX == 0x4003504f && ulHx > 0x333f9de6 ))){	// check qarea
		dfXsq = 1.0 / dfX;
		dfXsq = dfXsq*dfXsq;
		dfRet=K4*dfXsq+K3;
		dfRet=dfRet*dfXsq+K2;
		dfRet=dfRet*dfXsq+K1;
		dfRet=dfRet*dfXsq+K0;
		dfTemp=dfXsq+S3;
		dfTemp=dfTemp*dfXsq+S2;
		dfTemp=dfTemp*dfXsq+S1;
		dfTemp=dfTemp*dfXsq+S0;
		dfRet = HLFP-dfRet/dfTemp*(1.0/dfX);
		return dfRet;
	}else{
		dfXsq=(dfX-1.0)/(dfX+1.0);
		dfXsq=dfXsq*dfXsq;
		dfRet=K4*dfXsq+K3;
		dfRet=dfRet*dfXsq+K2;
		dfRet=dfRet*dfXsq+K1;
		dfRet=dfRet*dfXsq+K0;
		dfTemp=dfXsq+S3;
		dfTemp=dfTemp*dfXsq+S2;
		dfTemp=dfTemp*dfXsq+S1;
		dfTemp=dfTemp*dfXsq+S0;
		dfRet = QUAP+dfRet/dfTemp*((dfX - 1.0) / (dfX + 1.0));
		return dfRet;
	}
}